spatial_point_pattern

Load Packages

packages = c('rgdal', 'sf', 'tmap', 'tidyverse', 'sp', 'rgeos','maptools', 'raster', 'spatstat', 'tmaptools', 'spdep', 'OpenStreetMap')
for (p in packages){
  if (!require(p, character.only = T)){
    install.packages(p)
  }
  library(p, character.only = T)
}

## Loading required package: rgdal

## Loading required package: sp

## rgdal: version: 1.5-16, (SVN revision 1050)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: C:/Users/keith/Documents/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: C:/Users/keith/Documents/R/win-library/4.0/rgdal/proj
## Linking to sp version:1.4-2
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.

## Loading required package: sf

## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

## Loading required package: tmap

## Loading required package: tidyverse

## -- Attaching packages -------------------------------------- tidyverse 1.3.0 --

## v ggplot2 3.3.2     v purrr   0.3.4
## v tibble  3.0.3     v dplyr   1.0.1
## v tidyr   1.1.1     v stringr 1.4.0
## v readr   1.3.1     v forcats 0.5.0

## -- Conflicts ----------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

## Loading required package: rgeos

## rgeos version: 0.5-3, (SVN revision 634)
##  GEOS runtime version: 3.8.0-CAPI-1.13.1 
##  Linking to sp version: 1.4-2 
##  Polygon checking: TRUE

## Loading required package: maptools

## Checking rgeos availability: TRUE

## Loading required package: raster

## 
## Attaching package: 'raster'

## The following object is masked from 'package:dplyr':
## 
##     select

## The following object is masked from 'package:tidyr':
## 
##     extract

## Loading required package: spatstat

## Loading required package: spatstat.data

## Loading required package: nlme

## 
## Attaching package: 'nlme'

## The following object is masked from 'package:raster':
## 
##     getData

## The following object is masked from 'package:dplyr':
## 
##     collapse

## Loading required package: rpart

## Registered S3 method overwritten by 'spatstat':
##   method     from
##   print.boxx cli

## 
## spatstat 1.64-1       (nickname: 'Help you I can, yes!') 
## For an introduction to spatstat, type 'beginner'

## 
## Note: spatstat version 1.64-1 is out of date by more than 5 months; we recommend upgrading to the latest version.

## 
## Attaching package: 'spatstat'

## The following objects are masked from 'package:raster':
## 
##     area, rotate, shift

## Loading required package: tmaptools

## Loading required package: spdep

## Loading required package: spData

## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`

## Loading required package: OpenStreetMap

popdata <- read_csv('data/planning-area-subzone-age-group-sex-and-type-of-dwelling-june-2011-2019.csv')

## Parsed with column specification:
## cols(
##   planning_area = col_character(),
##   subzone = col_character(),
##   age_group = col_character(),
##   sex = col_character(),
##   type_of_dwelling = col_character(),
##   resident_count = col_double(),
##   year = col_double()
## )

popdata2019 <- popdata %>%
  filter(year == 2019) %>%
  filter(age_group == "65_to_69" | age_group == "70_to_74" | age_group == '75_to_79' | age_group == '80_to_84' | age_group == '85_to_89' | age_group == '90_and_over') %>%
  group_by(planning_area, subzone, age_group) %>%
  summarise(elderly_count = sum(resident_count)) %>%
  ungroup() %>%
  spread(age_group, elderly_count) %>%
  mutate(`elderly_count` = `65_to_69` + `70_to_74` + `75_to_79` + `80_to_84` + `85_to_89` + `90_and_over`)

## `summarise()` regrouping output by 'planning_area', 'subzone' (override with `.groups` argument)

popdata2019 <- mutate_at(popdata2019, .vars = c("subzone", "planning_area"), .funs=toupper)

Eldercare

eldercare_sf <- st_read(dsn='data', layer='ELDERCARE')

## Reading layer `ELDERCARE' from data source `C:\Users\keith\Desktop\IS415 Geospatial Analytics\SPP\SPP\data' using driver `ESRI Shapefile'
## Simple feature collection with 133 features and 18 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 14481.92 ymin: 28218.43 xmax: 41665.14 ymax: 46804.9
## projected CRS:  SVY21

eldercare_sf <- eldercare_sf %>%
  mutate(label = "Eldercare centres") 

Assign EPSG

eldercare_sf <- st_transform(eldercare_sf, 3414)

st_crs(eldercare_sf)

## Coordinate Reference System:
##   User input: EPSG:3414 
##   wkt:
## PROJCRS["SVY21 / Singapore TM",
##     BASEGEOGCRS["SVY21",
##         DATUM["SVY21",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4757]],
##     CONVERSION["Singapore Transverse Mercator",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",1.36666666666667,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",103.833333333333,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",1,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",28001.642,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",38744.572,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["northing (N)",north,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["easting (E)",east,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["Singapore"],
##         BBOX[1.13,103.59,1.47,104.07]],
##     ID["EPSG",3414]]

Convert sf to sp

eldercare_sp <- as_Spatial(eldercare_sf)

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on WGS84 ellipsoid in CRS definition,
##  but +towgs84= values preserved

eldercare_spatialpoint <- as(eldercare_sp, "SpatialPoints")

Silver Infocomm

infocomm_sf <- st_read(dsn='data', layer='SILVERINFOCOMM')

## Reading layer `SILVERINFOCOMM' from data source `C:\Users\keith\Desktop\IS415 Geospatial Analytics\SPP\SPP\data' using driver `ESRI Shapefile'
## Simple feature collection with 101 features and 18 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 12024.87 ymin: 28385.78 xmax: 41300.53 ymax: 47763.05
## projected CRS:  SVY21

infocomm_sf <- infocomm_sf %>%
  mutate(label = "Silver Infocomm Junc")

Assign EPSG

infocomm_sf <- st_transform(infocomm_sf, 3414)

st_crs(infocomm_sf)

## Coordinate Reference System:
##   User input: EPSG:3414 
##   wkt:
## PROJCRS["SVY21 / Singapore TM",
##     BASEGEOGCRS["SVY21",
##         DATUM["SVY21",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4757]],
##     CONVERSION["Singapore Transverse Mercator",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",1.36666666666667,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",103.833333333333,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",1,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",28001.642,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",38744.572,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["northing (N)",north,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["easting (E)",east,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["Singapore"],
##         BBOX[1.13,103.59,1.47,104.07]],
##     ID["EPSG",3414]]

Convert sf to sp

infocomm_sp <- as_Spatial(infocomm_sf)

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on WGS84 ellipsoid in CRS definition,
##  but +towgs84= values preserved

infocomm_spatialpoint <- as(infocomm_sp, "SpatialPoints")

CHAS Clinics

chas_sf <- st_read(dsn='data/chas-clinics-kml.kml')

## Reading layer `MOH_CHAS_CLINICS' from data source `C:\Users\keith\Desktop\IS415 Geospatial Analytics\SPP\SPP\data\chas-clinics-kml.kml' using driver `KML'
## Simple feature collection with 1167 features and 2 fields
## geometry type:  POINT
## dimension:      XYZ
## bbox:           xmin: 103.5818 ymin: 1.016264 xmax: 103.9903 ymax: 1.456037
## z_range:        zmin: 0 zmax: 0
## geographic CRS: WGS 84

chas_sf <- chas_sf %>%
  mutate(label = "Chas Clinics") %>%
  mutate(capacity = 1)

Assign EPSG

chas_sf <- st_transform(chas_sf, 3414)

st_crs(chas_sf)

## Coordinate Reference System:
##   User input: EPSG:3414 
##   wkt:
## PROJCRS["SVY21 / Singapore TM",
##     BASEGEOGCRS["SVY21",
##         DATUM["SVY21",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4757]],
##     CONVERSION["Singapore Transverse Mercator",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",1.36666666666667,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",103.833333333333,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",1,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",28001.642,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",38744.572,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["northing (N)",north,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["easting (E)",east,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["Singapore"],
##         BBOX[1.13,103.59,1.47,104.07]],
##     ID["EPSG",3414]]

Convert sf to sp

chas_sp <- as_Spatial(chas_sf)

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on WGS84 ellipsoid in CRS definition,
##  but +towgs84= values preserved

chas_spatialpoint <- as(chas_sp, "SpatialPoints")

MP_14 Subzone Data

mpsz_sf <- st_read(dsn='data', layer='MP14_SUBZONE_WEB_PL')

## Reading layer `MP14_SUBZONE_WEB_PL' from data source `C:\Users\keith\Desktop\IS415 Geospatial Analytics\SPP\SPP\data' using driver `ESRI Shapefile'
## Simple feature collection with 323 features and 15 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
## projected CRS:  SVY21

mpsz_sf <- mpsz_sf %>%
  dplyr::select(SUBZONE_N, PLN_AREA_N, REGION_N, X_ADDR, Y_ADDR, SHAPE_Leng, SHAPE_Area, geometry)

Assign EPSG

mpsz_sf <- st_transform(mpsz_sf, 3414)

st_crs(mpsz_sf)

## Coordinate Reference System:
##   User input: EPSG:3414 
##   wkt:
## PROJCRS["SVY21 / Singapore TM",
##     BASEGEOGCRS["SVY21",
##         DATUM["SVY21",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4757]],
##     CONVERSION["Singapore Transverse Mercator",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",1.36666666666667,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",103.833333333333,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",1,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",28001.642,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",38744.572,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["northing (N)",north,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["easting (E)",east,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["Singapore"],
##         BBOX[1.13,103.59,1.47,104.07]],
##     ID["EPSG",3414]]

centroids <- st_centroid(mpsz_sf)

## Warning in st_centroid.sf(mpsz_sf): st_centroid assumes attributes are constant
## over geometries of x

Coastal outline

sg <- readOGR(dsn = "data", layer="CostalOutline")

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\keith\Desktop\IS415 Geospatial Analytics\SPP\SPP\data", layer: "CostalOutline"
## with 60 features
## It has 4 fields

sg_spatialpoint <- as(sg, "SpatialPolygons")

Joining popdata2019 and mpsz_sf

mpsz_demand <- left_join(mpsz_sf, popdata2019, by=c('PLN_AREA_N' = 'planning_area', 'SUBZONE_N' = 'subzone'))
summary(mpsz_demand)

##   SUBZONE_N          PLN_AREA_N          REGION_N             X_ADDR     
##  Length:323         Length:323         Length:323         Min.   : 5093  
##  Class :character   Class :character   Class :character   1st Qu.:21864  
##  Mode  :character   Mode  :character   Mode  :character   Median :28465  
##                                                           Mean   :27257  
##                                                           3rd Qu.:31674  
##                                                           Max.   :50425  
##      Y_ADDR        SHAPE_Leng        SHAPE_Area          65_to_69     
##  Min.   :19579   Min.   :  871.5   Min.   :   39438   Min.   :   0.0  
##  1st Qu.:31776   1st Qu.: 3709.6   1st Qu.:  628261   1st Qu.:   0.0  
##  Median :35113   Median : 5211.9   Median : 1229894   Median : 270.0  
##  Mean   :36106   Mean   : 6524.4   Mean   : 2420882   Mean   : 686.7  
##  3rd Qu.:39869   3rd Qu.: 6942.6   3rd Qu.: 2106483   3rd Qu.:1060.0  
##  Max.   :49553   Max.   :68083.9   Max.   :69748299   Max.   :8120.0  
##     70_to_74         75_to_79         80_to_84         85_to_89  
##  Min.   :   0.0   Min.   :   0.0   Min.   :   0.0   Min.   :  0  
##  1st Qu.:   0.0   1st Qu.:   0.0   1st Qu.:   0.0   1st Qu.:  0  
##  Median : 180.0   Median : 110.0   Median :  70.0   Median : 40  
##  Mean   : 466.4   Mean   : 294.6   Mean   : 193.4   Mean   :105  
##  3rd Qu.: 760.0   3rd Qu.: 485.0   3rd Qu.: 310.0   3rd Qu.:170  
##  Max.   :4800.0   Max.   :2680.0   Max.   :1760.0   Max.   :970  
##   90_and_over     elderly_count            geometry  
##  Min.   :  0.00   Min.   :    0   MULTIPOLYGON :323  
##  1st Qu.:  0.00   1st Qu.:    0   epsg:3414    :  0  
##  Median : 30.00   Median :  730   +proj=tmer...:  0  
##  Mean   : 59.81   Mean   : 1806                      
##  3rd Qu.: 90.00   3rd Qu.: 3000                      
##  Max.   :530.00   Max.   :18860

mpsz_demand_sp <- as_Spatial(mpsz_demand)

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on WGS84 ellipsoid in CRS definition,
##  but +towgs84= values preserved

mpsz_demand_spatialpoint <- as(mpsz_demand_sp, "SpatialPolygons")

Calculating Chas Clinics count

mpsz_demand$`chas_count` <- lengths(st_intersects(mpsz_sf,chas_sf))

Calculating eldercare count

mpsz_demand$`eldercare_count` <- lengths(st_intersects(mpsz_sf,eldercare_sf))

Calculating infocomm count

mpsz_demand$`infocomm_count` <- lengths(st_intersects(mpsz_sf,infocomm_sf))

Convert sf to sp

mpsz_sp <- as_Spatial(mpsz_sf)

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on WGS84 ellipsoid in CRS definition,
##  but +towgs84= values preserved

mpsz_spatialpoint <- as(mpsz_sp, "SpatialPolygons")

sg_osm <- read_osm(mpsz_spatialpoint, ext=1.3)

Class interval for boxmap

boxbreaks <- function(v,mult=1.5) {
  qv <- unname(quantile(v))
  iqr <- qv[4] - qv[2]
  upfence <- qv[4] + mult * iqr
  lofence <- qv[2] - mult * iqr
  # initialize break points vector. might need more length for bigger data
  bb <- vector(mode="numeric",length=7)
  # logic for lower and upper fences
  if (lofence < qv[1]) { # no lower outliers
    bb[1] <- lofence
    bb[2] <- floor(qv[1])
  } else {
    bb[2] <- lofence
    bb[1] <- qv[1]
  }
  if (upfence > qv[5]) { # no upper outliers
    bb[7] <- upfence
    bb[6] <- ceiling(qv[5])
  } else {
    bb[6] <- upfence
    bb[7] <- qv[5]
  }
  bb[3:5] <- qv[2:4]
  return(bb)
}

get.var <- function(vname,df) {
  v <- df[vname] %>% st_set_geometry(NULL)
  v <- unname(v[,1])
  return(v)
}

boxmap <- function(vnam,df,legtitle=NA,mtitle="Box Map",mult=1.5){
  var <- get.var(vnam,df)
  bb <- boxbreaks(var)
  tm_shape(df) +
    tm_fill(vnam,title=legtitle,breaks=bb,palette="-RdBu",
            labels = c("lower outlier", "< 25%", "25% - 50%", "50% - 75%","> 75%", "upper outlier")) +
  tm_borders() +
  tm_layout(title = mtitle, title.position = c("right","bottom"))
}

boxmap("elderly_count", mpsz_demand, mtitle="Elderly Distribution")

Demand and supply of eldercare facilities

tmap_mode("view")

## tmap mode set to interactive viewing

tm_shape(mpsz_demand) + 
  tm_fill(c("elderly_count"),
          title = "Elderly Population",
          breaks = c(0, 2500 ,5000, 7500, 10000, 12500, 15000, 17500, 20000),
          palette = "Blues",
          popup.vars=c("Subzone Name"="SUBZONE_N", 
                       "Planning Area Name"="PLN_AREA_N", 
                       "Elderly Population"="elderly_count"),
          showNA = FALSE) +
  tm_borders(lwd = 0.1,  alpha = 1) +
  tm_shape(chas_sf) + tm_dots("red") +
  tm_shape(eldercare_sf) + tm_dots("blue")+
  tm_shape(infocomm_sf) + tm_dots("green") 

Extracting study areas(Tampines, Ang Mo Kio, Serangoon, Jurong West) in spatialpoint dataframe

tp <- mpsz_sp[mpsz_sp@data$PLN_AREA_N == "TAMPINES",]
amk <- mpsz_sp[mpsz_sp@data$PLN_AREA_N == "ANG MO KIO",]
bs <- mpsz_sp[mpsz_sp@data$PLN_AREA_N == "BISHAN",]
jw <- mpsz_sp[mpsz_sp@data$PLN_AREA_N == "JURONG WEST",]

Convert them into generic spatialpolygons objects

tp_spatialpoint <- as(tp, "SpatialPolygons")
amk_spatialpoint <- as(amk, "SpatialPolygons")
bs_spatialpoint <- as(bs, "SpatialPolygons")
jw_spatialpoint <- as(jw, "SpatialPolygons")

creating own objects based on subzones

tp_owin <- as(tp, "owin")
amk_owin <- as(amk, "owin")
bs_owin <- as(bs, "owin")
jw_owin <- as(jw, "owin")

Converting spatialpoints to spatstat ppp

chas_ppp <- as(chas_spatialpoint, "ppp")
infocomm_ppp <- as(infocomm_spatialpoint, "ppp")
eldercare_ppp <- as(eldercare_spatialpoint, "ppp")

using jittering to better visualize overlapped data

chas_ppp_jit <- rjitter(chas_ppp, retry=TRUE, nsim=1, drop=TRUE)
infocomm_ppp_jit <- rjitter(infocomm_ppp, retry=TRUE, nsim=1, drop=TRUE)
eldercare_ppp_jit <- rjitter(eldercare_ppp, retry=TRUE, nsim=1, drop=TRUE)

Combining Chas, Infocomm and Eldercare

ppp_chas_tp <- chas_ppp_jit[tp_owin]
ppp_chas_amk <- chas_ppp_jit[amk_owin]
ppp_chas_bs <- chas_ppp_jit[bs_owin]
ppp_chas_jw <- chas_ppp_jit[jw_owin]

ppp_infocomm_tp <- infocomm_ppp_jit[tp_owin]
ppp_infocomm_amk <- infocomm_ppp_jit[amk_owin]
ppp_infocomm_bs <- infocomm_ppp_jit[bs_owin]
ppp_infocomm_jw <- infocomm_ppp_jit[jw_owin]

ppp_eldercare_tp <- eldercare_ppp_jit[tp_owin]
ppp_eldercare_amk <- eldercare_ppp_jit[amk_owin]
ppp_eldercare_bs <- eldercare_ppp_jit[bs_owin]
ppp_eldercare_jw <- eldercare_ppp_jit[jw_owin]

Convert data to KM by rescaling

ppp_chas_tp.km = rescale(ppp_chas_tp, 1000, "km")
ppp_chas_amk.km = rescale(ppp_chas_amk, 1000, "km")
ppp_chas_bs.km = rescale(ppp_chas_bs, 1000, "km")
ppp_chas_jw.km = rescale(ppp_chas_jw, 1000, "km")

ppp_infocomm_tp.km = rescale(ppp_infocomm_tp, 1000, "km")
ppp_infocomm_amk.km = rescale(ppp_infocomm_amk, 1000, "km")
ppp_infocomm_bs.km = rescale(ppp_infocomm_bs, 1000, "km")
ppp_infocomm_jw.km = rescale(ppp_infocomm_jw, 1000, "km")

ppp_eldercare_tp.km = rescale(ppp_eldercare_tp, 1000, "km")
ppp_eldercare_amk.km = rescale(ppp_eldercare_amk, 1000, "km")
ppp_eldercare_bs.km = rescale(ppp_eldercare_bs, 1000, "km")
ppp_eldercare_jw.km = rescale(ppp_eldercare_jw, 1000, "km")

create owin objects that is required by spatstat.

sg_owin <- as(mpsz_sp, "owin")
sg_owin

## window: polygonal boundary
## enclosing rectangle: [2667.54, 56396.44] x [15748.72, 50256.33] units

Extracting all facilities within the boundaries of the study area

ppp_chas_mpsz = chas_ppp_jit[sg_owin]
ppp_infocomm_mpsz = infocomm_ppp_jit[sg_owin]
ppp_eldercare_mpsz = eldercare_ppp_jit[sg_owin]

Rescaling in KM

ppp_chas_mpsz.km = rescale(ppp_chas_mpsz, 1000, "km")
ppp_infocomm_mpsz.km = rescale(ppp_infocomm_mpsz, 1000, "km")
ppp_eldercare_mpsz.km = rescale(ppp_eldercare_mpsz, 1000, "km")

Spatial point pattern analysis for chas clinics

In order to confirm the spatial patterns for chas clinics, we will have to conduct a hypothesis test for all subzones and chas clinic locations.

H0 = The distribution of chas clinics are randomly distributed. H1= The distribution of chas clinics are not randomly distributed.

The null hypothesis will be rejected if p-value is smaller than alpha < 0.001i

converting spatialdata into ppp object format

performing CSR test

K_chas_tp = Lest(ppp_chas_tp, correction = "Ripley")
K_chas_amk = Lest(ppp_chas_amk, correction = "Ripley")
K_chas_bs = Lest(ppp_chas_bs, correction = "Ripley")
K_chas_jw = Lest(ppp_chas_jw, correction = "Ripley")

K_infocomm_tp = Lest(ppp_infocomm_tp, correction = "Ripley")
K_infocomm_amk = Lest(ppp_infocomm_amk, correction = "Ripley")
K_infocomm_bs = Lest(ppp_infocomm_bs, correction = "Ripley")
K_infocomm_jw = Lest(ppp_infocomm_jw, correction = "Ripley")

K_eldercare_tp = Lest(ppp_eldercare_tp, correction = "Ripley")
K_eldercare_amk = Lest(ppp_eldercare_tp, correction = "Ripley")
K_eldercare_bs = Lest(ppp_eldercare_tp, correction = "Ripley")
K_eldercare_jw = Lest(ppp_eldercare_tp, correction = "Ripley")

K_chas_tp.csr <- envelope(ppp_chas_tp, Kest, nsim = 99, rank = 1, glocal=TRUE)

## Generating 99 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.

K_chas_amk.csr <- envelope(ppp_chas_amk, Kest, nsim = 99, rank = 1, glocal=TRUE)

## Generating 99 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.

K_chas_bs.csr <- envelope(ppp_chas_bs, Kest, nsim = 99, rank = 1, glocal=TRUE)

## Generating 99 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.

K_chas_jw.csr <- envelope(ppp_chas_jw, Kest, nsim = 99, rank = 1, glocal=TRUE)

## Generating 99 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.

K_infocomm_tp.csr <- envelope(ppp_infocomm_tp, Kest, nsim = 99, rank = 1, glocal=TRUE)

## Generating 99 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.

K_infocomm_amk.csr <- envelope(ppp_infocomm_amk, Kest, nsim = 99, rank = 1, glocal=TRUE)

## Generating 99 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.

K_infocomm_bs.csr <- envelope(ppp_infocomm_bs, Kest, nsim = 99, rank = 1, glocal=TRUE)

## Generating 99 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.

K_infocomm_jw.csr <- envelope(ppp_infocomm_jw, Kest, nsim = 99, rank = 1, glocal=TRUE)

## Generating 99 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.

K_eldercare_tp.csr <- envelope(ppp_eldercare_tp, Kest, nsim = 99, rank = 1, glocal=TRUE)

## Generating 99 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.

K_eldercare_amk.csr <- envelope(ppp_eldercare_tp, Kest, nsim = 99, rank = 1, glocal=TRUE)

## Generating 99 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.

K_eldercare_bs.csr <- envelope(ppp_eldercare_tp, Kest, nsim = 99, rank = 1, glocal=TRUE)

## Generating 99 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.

K_eldercare_jw.csr <- envelope(ppp_eldercare_tp, Kest, nsim = 99, rank = 1, glocal=TRUE)

## Generating 99 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.

par(mfrow=c(4,3))
par(mar = rep(3, 4))

plot(K_chas_tp, . - r ~ r, 
     xlab="d", ylab="L(d)-r", 
     main="Chas Clinics in Tampines")
plot(K_infocomm_tp.csr, . - r ~ r, 
     xlab="d", ylab="L(d)-r", 
     main="Silver infocomm centres in Tampines")
plot(K_eldercare_tp.csr, . - r ~ r, 
     xlab="d", ylab="L(d)-r", 
     main="Eldercare services in Tampines")

plot(K_chas_amk.csr, . - r ~ r, 
     xlab="d", ylab="L(d)-r", 
     main="Chas Clinics in Ang Mo Kio")
plot(K_infocomm_amk.csr, . - r ~ r, 
     xlab="d", ylab="L(d)-r", 
     main="Silver infocomm centres in Ang Mo Kio")
plot(K_eldercare_amk.csr, . - r ~ r, 
     xlab="d", ylab="L(d)-r", 
     main="Eldercare services in Ang Mo Kio")

plot(K_chas_bs.csr, . - r ~ r, 
     xlab="d", ylab="L(d)-r", 
     main="Chas Clinics in Bishan")
plot(K_infocomm_bs.csr, . - r ~ r, 
     xlab="d", ylab="L(d)-r", 
     main="Silver infocomm centres in Bishan")
plot(K_eldercare_bs.csr, . - r ~ r, 
     xlab="d", ylab="L(d)-r", 
     main="Eldercare services in Bishan")

plot(K_chas_jw.csr, . - r ~ r, 
     xlab="d", ylab="L(d)-r", 
     main="Chas Clinics in Jurong West")
plot(K_infocomm_jw.csr, . - r ~ r, 
     xlab="d", ylab="L(d)-r", 
     main="Silver infocomm centres in Jurong West")
plot(K_eldercare_jw.csr, . - r ~ r, 
     xlab="d", ylab="L(d)-r", 
     main="Eldercare services in Jurong West")

##Kernal density Estimation with OSM

Create a binding box for selected subzones

tp_bb <- st_bbox(mpsz_sf %>% filter(PLN_AREA_N == "TAMPINES"))
amk_bb <- st_bbox(mpsz_sf %>% filter(PLN_AREA_N == "ANG MO KIO"))
bs_bb <- st_bbox(mpsz_sf %>% filter(PLN_AREA_N == "BISHAN"))
jw_bb <- st_bbox(mpsz_sf %>% filter(PLN_AREA_N == "JURONG WEST"))

Getting osm for the indivisual planning areas

tp_osm <- read_osm(tp_bb, ext=1.1)
amk_osm <- read_osm(amk_bb, ext=1.1)
bs_osm <- read_osm(bs_bb, ext=1.1)
jw_osm <- read_osm(jw_bb, ext=1.1)

getPlnAreaKernelDensityMap <- function(osm, pln, ppp, ppp_str) {
  ppp.km <- rescale(ppp, 1000, "km")
  
  kde <- density(ppp.km, sigma=0.20, edge=TRUE, kernel="gaussian")
  
  gridded_kde_bw <- as.SpatialGridDataFrame.im(kde)
  
  kde_bw_raster <- raster(gridded_kde_bw)
  
  projection(kde_bw_raster) <- crs("+init=EPSG:3414 +datum=WGS84 +units=km")
  
  # Plot kernel density map on openstreetmap
  tm_shape(osm)+ 
    tm_layout(legend.outside = TRUE, title=ppp_str)+
    tm_rgb()+
  tm_shape(pln)+
    tm_borders(col = "darkblue", lwd = 2, lty="longdash")+
  tm_shape(kde_bw_raster) + 
    tm_raster("v", alpha=0.5,  
          palette = "YlOrRd")
  
}

kde_chas_tp <- getPlnAreaKernelDensityMap(tp_osm, tp, ppp_chas_tp, "Chas Clinics in Tampines") 
kde_infocomm_tp <- getPlnAreaKernelDensityMap(tp_osm, tp, ppp_infocomm_tp, "Infocomm Centres in Tampines")
kde_eldercare_tp <- getPlnAreaKernelDensityMap(tp_osm, tp, ppp_eldercare_tp, "Eldercare Services in Tampines")

kde_chas_amk <- getPlnAreaKernelDensityMap(amk_osm, amk, ppp_chas_amk, "Chas Clinics in Ang Mo Kio") 
kde_infocomm_amk <- getPlnAreaKernelDensityMap(amk_osm, amk, ppp_infocomm_amk, "Infocomm Centres in Ang Mo Kio")
kde_eldercare_amk <- getPlnAreaKernelDensityMap(amk_osm, amk, ppp_eldercare_amk, "Eldercare Services in Ang Mo Kio")

kde_chas_bs <- getPlnAreaKernelDensityMap(bs_osm, bs, ppp_chas_bs, "Chas Clinics in Bishan") 
kde_infocomm_bs <- getPlnAreaKernelDensityMap(bs_osm, bs, ppp_infocomm_bs, "Infocomm Centres in Bishan")
kde_eldercare_bs <- getPlnAreaKernelDensityMap(bs_osm, bs, ppp_eldercare_bs, "Eldercare Services in Bishan")

kde_chas_jw <- getPlnAreaKernelDensityMap(jw_osm, jw, ppp_chas_jw, "Chas Clinics in Jurong West") 
kde_infocomm_jw <- getPlnAreaKernelDensityMap(jw_osm, jw, ppp_infocomm_jw, "Infocomm Centres in Jurong West")
kde_eldercare_jw <- getPlnAreaKernelDensityMap(jw_osm, jw, ppp_eldercare_jw, "Eldercare Services in Jurong West")

tmap_arrange(kde_chas_tp, kde_infocomm_tp, kde_eldercare_tp, asp=1, ncol=3)

## OpenStreetMapData read with read_osm is static, so not usable in view mode. Please use tm_basemap or tm_tiles, with the provider name set to "OpenStreetMap.Mapnik"
## OpenStreetMapData read with read_osm is static, so not usable in view mode. Please use tm_basemap or tm_tiles, with the provider name set to "OpenStreetMap.Mapnik"
## OpenStreetMapData read with read_osm is static, so not usable in view mode. Please use tm_basemap or tm_tiles, with the provider name set to "OpenStreetMap.Mapnik"

tmap_arrange(kde_chas_amk, kde_infocomm_amk, kde_eldercare_amk, asp=1, ncol=3)

## OpenStreetMapData read with read_osm is static, so not usable in view mode. Please use tm_basemap or tm_tiles, with the provider name set to "OpenStreetMap.Mapnik"
## OpenStreetMapData read with read_osm is static, so not usable in view mode. Please use tm_basemap or tm_tiles, with the provider name set to "OpenStreetMap.Mapnik"
## OpenStreetMapData read with read_osm is static, so not usable in view mode. Please use tm_basemap or tm_tiles, with the provider name set to "OpenStreetMap.Mapnik"

tmap_arrange(kde_chas_bs, kde_infocomm_bs, kde_eldercare_bs, asp=1, ncol=3)

## OpenStreetMapData read with read_osm is static, so not usable in view mode. Please use tm_basemap or tm_tiles, with the provider name set to "OpenStreetMap.Mapnik"
## OpenStreetMapData read with read_osm is static, so not usable in view mode. Please use tm_basemap or tm_tiles, with the provider name set to "OpenStreetMap.Mapnik"
## OpenStreetMapData read with read_osm is static, so not usable in view mode. Please use tm_basemap or tm_tiles, with the provider name set to "OpenStreetMap.Mapnik"

tmap_arrange(kde_chas_jw, kde_infocomm_jw, kde_eldercare_jw, asp=1, ncol=3)

## OpenStreetMapData read with read_osm is static, so not usable in view mode. Please use tm_basemap or tm_tiles, with the provider name set to "OpenStreetMap.Mapnik"
## OpenStreetMapData read with read_osm is static, so not usable in view mode. Please use tm_basemap or tm_tiles, with the provider name set to "OpenStreetMap.Mapnik"
## OpenStreetMapData read with read_osm is static, so not usable in view mode. Please use tm_basemap or tm_tiles, with the provider name set to "OpenStreetMap.Mapnik"